Cone beam local tomography

ABSTRACT

Methods, systems and processes for providing efficient image reconstruction using local cone beam tomography which provide a reduced level of artifacts without suppressing the strength of the useful features; and in a dynamic case provide reconstruction of objects that are undergoing a change during the scan. An embodiment provides a method of reconstructing an image from cone beam data provided by at least one detector. The method includes collecting CB projection data of an object, storing the CB projection data in a memory; and reconstructing the image from the local CB projection data. In the reconstructing step, a combination of derivatives of the CB projection data that will result in suppressing the artifacts are found. The combination of derivatives includes collecting cone beam data that represents a collection of integrals that represent the object.

This Patent Application claims the benefit of priority to U.S.Provisional Patent Application No. 60/754,236 filed on Dec. 28, 2005 andwas funded in part by NSF grant DMS-0104033 and the Alexander vonHumboldt Foundation.

FIELD OF THE INVENTION

This invention relates to local tomography and, in particular, tomethods, systems, apparatus and devices for cone beam local tomography.

BACKGROUND AND PRIOR ART

In computed tomography (CT) the goal is to reconstruct the distributionof the x-ray attenuation coefficient ƒ inside the object being scanned.Local tomography (LT) computes not ƒ, but Bƒ, where B is some operatorthat enhances singularities of ƒ. In two dimensions (2D), B is anelliptic pseudo-differential operator (PDO) of order one as described byA. Katsevich, “Local Tomography for the Generalized Random Transform”,SIAM Journal on Applied Mathematics, Vol. 57, no. 4 (1997) pp.1128-1162, Kuchment et al., “On local Tomography”, Inverse Problem 11(1995), pp. 571-589, A. Ramm and A. Katsevich, “The Random Transform andLocal Tomography”, CRC Press, Boca Raton, Fla., 1996 and A. G. Ramm,“Necessary and Sufficient Conditions for a PDO to be a local tomographyoperator”, Comptes Rend Acad. Sci., Paris 332 (1996) pp. 613-618. In thecone beam setting (three dimensions) a LT function is denoted by g_(Λ).The corresponding operator B: ƒ→g_(Λ) is much more complicated than in2D. It preserves the so-called visible (or, useful) singularities andcreates non-local artifacts. Unfortunately, the strength of theseartifacts is the same as that of the useful singularities of g_(Λ) asdescribed in A. Katsevich “Cone Beam Local Tomography”, SIAM Journal onApplied Mathematics (1999), pp. 2224-2246 and D. Finch et al,“Microlocal analysis of the X-ray transform with sources on a curve”,Inside out: Inverse Problems and Applications, Cambridge Univ. Press,(2003) pp. 193-218.

SUMMARY OF THE INVENTION

A primary objective of the invention is to provide a new method, system,apparatus and device for improved cone beam local tomography.

A secondary objective of the invention is to provide new methods,systems, apparatus and devices for cone beam local tomography thatproduces artifacts that are one order smoother in the scale of Sobolevspaces.

A third objective of the invention is to provide new methods, systems,apparatus and devices for cone beam local tomography using theflexibility of LT, its relative stability with respect toinconsistencies in the data and the ability to accurately reconstructedges inside objects, which provides important information complementingwell-established inversion techniques.

A fourth objective of the invention is to provide methods, systems,apparatus, and devices for reconstructing an image from local cone beamprojection data with suppressed artifacts.

A fifth objective of the invention is to provide methods, systems,apparatus, and devices for reconstructing an image from local cone beamprojection data with suppressed artifacts without suppressing the usefulfeatures of the image.

A sixth objective of the invention is to provide methods, systems,apparatus, and devices for determining a shift between a location of thereconstruction point shown on the reconstructed image and an actuallocation of the point to determine a location error of the moving objectto correct the location of the useful features of the actual image.

A first embodiment provides a method of reconstructing an image fromcone beam data provided by at least one detector. First, cone beamprojection data of an object is collected and stored in memory. Theimage is reconstructed from the local cone beam projection data, whereinartifacts of the image are suppressed without suppressing the strengthof the useful features. During the reconstructing step, a direction of aderivative of the cone beam projection data that will result insuppressing the artifacts is found, the direction being a tangent to thecurve, which represents the scan trajectory. Based on the direction, thecone beam projection data is differentiated along the direction of thetangent of the curve. The derivative results are processed, then backprojection is applied to the processed data to extract the usefulfeatures of the image with suppressed artifacts.

In an embodiment, a shift between a location of the reconstruction pointshown on the reconstructed image and an actual location of the point inthe object is determined to correct for the location of the usefulfeatures of the actual image.

Further objects and advantages of this invention will be apparent fromthe following detailed description of preferred embodiments which areillustrated schematically in the accompanying drawings.

BRIEF DESCRIPTION OF THE FIGURES

FIG. 1 shows an overview of the basic process steps according to theinvention.

FIG. 2 is flow diagram for the back-projection substeps corresponding tostep 40 of FIG. 1.

FIG. 3 a shows the new LT function g in the case of the clock phantom.

FIG. 3 b shows old LT function g_(Λ) in the case of the clock phantom.

FIG. 4 a shows the results of reconstructing the three-disk phantomaccording to the present invention.

FIG. 4 b shows the results of reconstructing the three-disk phantomaccording to the prior art.

FIG. 5 a shows the results for a phantom consisting of two identicalballs having radius 40 centered at (−80,0,0) and (80,0,0), respectively,according to the present invention.

FIG. 5 b shows the results for a phantom consisting of two identicalballs having radius 40 centered at (−80,0,0) and (80,0,0), respectively,according to the prior art.

FIG. 6 a shows reconstruction of an ellipsoid in the dynamic case, newLT function g according to the present invention wherein the ellipsoidmoves either up or down.

FIG. 6 b shows reconstruction of an ellipsoid in the dynamic case, newLT function g according to the present invention wherein the ellipsoidmoves either left or right

FIG. 7 a shows reconstruction of an ellipsoid in the dynamic case, newLT function g according to the present invention wherein the ellipsoidexpands and contracts.

FIG. 7 b shows reconstruction of an ellipsoid in the dynamic case, theprior art LT function g wherein the ellipsoid expands and contracts aswell as moves left-right.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

Before explaining the disclosed embodiments of the present invention indetail it is to be understood that the invention is not limited in itsapplication to the details of the particular arrangements shown sincethe invention is capable of other embodiments. Also, the terminologyused herein is for the purpose of description and not of limitation.

Local Tomography Function of the Present Invention:

-   -   C is a smooth curve defined by I        s→y(s)⊂R³. Here    -   s is a real parameter;    -   I is a compact interval;    -   y(s) denote the parametric equations of the curve; and    -   R³ is the three-dimensional Euclidean space.        Fix an open set V, whose closure V⊂R³\C. Let D_(ƒ)(s,α) denote        the cone beam transform of ƒ:

$\begin{matrix}{{{D_{f}\left( {s,\alpha} \right)} = {\int_{0}^{\infty}{{f\left( {{y(s)} + {\alpha\; t}} \right)}{\mathbb{d}t}}}},{s \in I},{\alpha \in S^{2}},} & (1)\end{matrix}$where ƒ∈C₀ ^(∞)(V). Here

-   -   α is a unit vector;    -   S² is the unit sphere in R³; and    -   ƒ is the function representing the distribution of the x-ray        attenuation coefficient inside the object being scanned        Fix φ∈C₀ ^(∞)(I×V) and define:

$\begin{matrix}{{{g(x)}\text{:}{= {{B\;{f(x)}\text{:} = {\int_{I}{{\varphi\left( {s,x} \right)}\frac{\partial^{2}}{\partial q^{2}}{D_{f}\left( {s,{\beta\left( {q,x} \right)}} \right)}}}}❘_{q = s}{\mathbb{d}s}}}},} & (2)\end{matrix}$where q is real parameter and

$\begin{matrix}{{\beta\left( {q,x} \right)}\text{:}{= \frac{x - {y(q)}}{{x - {y(q)}}}.}} & (3)\end{matrix}$g(x) defined by equation (2) above is one example of a new localtomography function of the present invention. Introduce the followingnotationΠ(x,ξ):={y∈R ³:ξ·(y−x)=0},WF _(v)(ƒ):={(x,ξ)∈WF(ƒ):Π(x,ξ) intersects C transversely},L ₊(s,z):={x∈R ³ :x=y(s)+t(z−y(s)),t>0}.  (4)Here

-   -   ξ is a non-zero vector in R³;    -   z is a point in R³;    -   WF(ƒ) is the wave-front of ƒ; and    -   WF_(v)(ƒ) denotes the visible singularities of ƒ. The set        WF_(v)(ƒ) and the term ‘visible singularity’ were introduced        in E. T. Quinto, “Singularies of the X-ray transform and limited        data tomography in R² and R³”, SIAM Journal of Mathematical        Analysis, vol. 24 (1993), pp. 1215-1225.        Analogously to A. Katsevich, “Cone beam local tomography”, SIAM        Journal on Applied Mathematics, vol. 59 (1999), pp. 2224-2246,        the following results can be established: the operator B defined        by equation (2) extends to a map E′(V)→E′(V), and        WF(g)⊂WF_(v)(ƒ)∪E(ƒ,C),        E(ƒ,C):={(x,ξ)∈T*V\0:ξ⊥{dot over (y)}(s),ξ⊥(x−y(s)),x∈L        ₊(s,z),(z,ξ)∈WF(ƒ),(s,x)∈suppφ}.  (5)        Here E′(V) is the space of distributions with compact support        inside V, and V is an open set not intersecting C.

An equivalent definition of the set E(ƒ,C) is as follows. If (z,ξ)∈WF(ƒ)and the plane Π(z,ξ) is tangent to C at some point y(s), then(x,ξ)∈E(ƒ,C) for all x∈L₊(s,z), (s,x)∈suppφ.

ξ⁽⁰⁾≠0 and x⁽⁰⁾∈R³\C are selected such that the plane Π(x⁽⁰⁾,ξ⁽⁰⁾)intersects C transversely. Then there exists a sufficiently small conicneighborhood U×Ω⊂T*V\0 of (x⁽⁰⁾,ξ⁽⁰⁾) such that B(x,ξ) is a classicalamplitude from the class S¹(U×Ω). Moreover, if Π(x⁽⁰⁾,ξ⁽⁰⁾) does notintersect C, then B(x,ξ) is from the class S^(−∞)(U×Ω) for someneighborhood U×Ω

(x⁽⁰⁾,ξ⁽⁰⁾).

The principal symbol of B is computed as follows. Denote

$\begin{matrix}{{m\text{:} = \inf\limits_{{({s,x})} \in {{supp}\;\varphi}}{{x - {y(s)}}}},{M\text{:} = \sup\limits_{{({s,x})} \in {{supp}\;\varphi}}{{x - {y(s)}}}},} & (6)\end{matrix}$and pick δ, 0<δ<m. Let w(t) be a function with the propertiesw(t)∈C ₀ ^(∞)([m−δ,M+δ]),w(t)=1,t∈[m,M].  (7)Representing ƒ in terms of its Fourier transform and using equation (7)we get from equations (1) and (2)

$\begin{matrix}{{{g(x)} = {{{\frac{1}{\left( {2\pi} \right)^{3}}{\int_{R^{3}}{{\overset{\sim}{f}(\xi)}{\int_{R^{2}}{{\varphi\left( {s,x} \right)} \times \frac{\partial^{2}}{\partial q^{2}}{w\left( {t{{x - {y(q)}}}} \right)}{{x - {y(q)}}}{\mathbb{e}}^{{- {\mathbb{i}\xi}} \cdot {({{y{(s)}} + {t{({x - {y{(q)}}})}}})}}}}}}}❘_{q = s}{{\mathbb{d}t}{\mathbb{d}s}{\mathbb{d}\xi}}} = {\frac{1}{\left( {2\pi} \right)^{3}}{\int_{R^{3}}{{\overset{\sim}{f}(\xi)}{B\left( {x,\xi} \right)}{\mathbb{e}}^{{- {\mathbb{i}\xi}} \cdot x}{\mathbb{d}\xi}}}}}},{where}} & (8) \\{{B\left( {x,\xi} \right)} = {- {\int_{R^{2}}{{{\varphi\left( {s,x} \right)}\left\lbrack {{{w\left( {t{{x - {y(s)}}}} \right)}{{x - {y(s)}}}\left( {t\;{\xi \cdot {\overset{.}{y}(s)}}} \right)^{2}} + {O(\xi)}} \right\rbrack} \times {\mathbb{e}}^{{{- {\mathbb{i}\xi}} \cdot {({x - {y{(s)}}})}}{({t - 1})}}{\mathbb{d}t}{{\mathbb{d}s}.}}}}} & (9)\end{matrix}$Pick (x,ξ)∈T*V\0 such that Π(x,ξ) intersects C, and the intersectionsare transversal. Asymptotics of B(x,σξ) as σ→∞ are found. Here σ is areal parameter. The stationary points of the phase in equation (9) aredetermined by solving(ξ·{dot over (y)}(s))(t−1)=0, ξ·(x−y(s))=0  (10)for s and t. Let s_(j)=s_(j)(x,ξ) be the points of intersections of theplane Π(x,ξ) and curve C. As is easily seen, the critical points are(s=s_(j),t=1). They are non-degenerate, because by assumption ξ·{dotover (y)}(s_(j))≠0. Now the stationary phase method immediately yields

$\begin{matrix}{{{B\left( {x,{\sigma\xi}} \right)} = {{{- 2}{\pi\sigma}{\sum\limits_{j}{{\varphi\left( {s_{j},x} \right)}{{x - {y\left( s_{j} \right)}}}{{\xi \cdot {\overset{.}{y}\left( s_{j} \right)}}}}}} + {O(1)}}},{\sigma->{\infty.}}} & (11)\end{matrix}$Analysis of Artifacts:

Pick any (z⁽⁰⁾,ξ⁽⁰⁾)∈WF(ƒ) such that the plane Π(z⁽⁰⁾,ξ⁽⁰⁾) is tangentto C at some y(s₀). Fix x⁽⁰⁾∈L₊(s₀,z⁽⁰⁾), x⁽⁰⁾∉{z⁽⁰⁾,y(s₀)} and(s₀,x⁽⁰⁾)∈suppφ. The problem is to find the behavior of g, defined byequation (2), in a neighborhood of x⁽⁰⁾. This behavior, of course,depends on the nature of the singularity of ƒ at z⁽⁰⁾. Assuming that ƒis a conormal distribution associated to a smooth hypersurface S, whichhas nonzero Gaussian curvature at z⁽⁰⁾: ƒ∈I_(comp) ^(m)(V,N*S) for somem. Here N*S is the conormal bundle of S. Let U×Ω be a sufficiently smallconic neighborhood of (z⁽⁰⁾,ξ⁽⁰⁾). In view of the above examples, it isassumed without loss of generality that {tilde over(ƒ)}(ξ)=A(ξ)e^(iH(ξ))+A₁(ξ) (cf. L. Hormander, The Analysis of LinearPartial Differential Operators, Volume IV, Springer Verlag, New York,1985, Proposition 25.1.3). Here A∈S^(m−3/4)(R³) and A(ξ)≡0 on R³\Ω,A₁∈S^(−∞)(R³), H(ξ)∈C₀ ^(∞)(R³\0) is real-valued and homogeneous ofdegree one, and N*S=(H′(ξ),ξ) on U×Ω.

The coordinate system is introduced with the origin at z⁽⁰⁾, such thatξ⁽⁰⁾ is along the x₃-axis, and y(s₀)=(c,0,0), where c=|y(s₀)−z⁽⁰⁾|. Then

$\begin{matrix}{{{x^{(0)} = \left( {b,0,0} \right)},{{{\overset{.}{y}}_{3}\left( s_{0} \right)} = 0},{{{{\overset{.}{y}}_{2}\left( s_{0} \right)} \neq 0};}}{{{H_{3}\left( \xi^{(0)} \right)} = {{H_{j\; 3}\left( \xi^{(0)} \right)} = 0}},{j = 1},2,{3;}}{{\begin{matrix}{H_{11}\left( \xi^{(0)} \right)} & {H_{12}\left( \xi^{(0)} \right)} \\{H_{21}\left( \xi^{(0)} \right)} & {H_{22}\left( \xi^{(0)} \right)}\end{matrix}} \neq 0.}} & (12)\end{matrix}$The inequality {dot over (y)}₂(s₀)≠0 is equivalent to the assumptionthat the line tangent to C at y(s₀) does not contain z⁽⁰⁾. Thisassumption is not very restrictive, because almost all scanningprotocols satisfy it. By assumption, x⁽⁰⁾∉{z⁽⁰⁾,y(s₀)}, so b∉{0,c}.Denote φ₁(s,x):=−φ(s,x)|x−y(s)|/(2π)³. Substituting the expression for{tilde over (ƒ)} into equation (8), using equation (9) produces

$\begin{matrix}{{{{g(x)} - {\int_{R^{2}}{{\varphi_{1}\left( {s,x} \right)}{\int_{\Omega}{{{A(\xi)}\left\lbrack {{{w\left( {t{{x - {y(s)}}}} \right)}\left( {t\;{\xi \cdot {\overset{.}{y}(s)}}} \right)^{2}} + {O(\xi)}} \right\rbrack} \times {\mathbb{e}}^{{\mathbb{i}}{({{H{(\xi)}} - {\xi \cdot {z{({s,t})}}}})}}{\mathbb{d}\xi}{\mathbb{d}t}{\mathbb{d}s}}}}}} \in {C_{0}^{\infty}(V)}},\mspace{20mu}{{{z\left( {s,t} \right)}\text{:}} = {{y(s)} + {t\left( {x - {y(s)}} \right)}}},} & (13)\end{matrix}$where the integral with respect to ξ is understood as an oscillatoryone. For simplicity of notation, the dependence of z(s,t) on x isomitted. Let I denote the integral with respect to ξ in equation (13).Changing variables ξ→({circumflex over (ξ)},λ), where ξ=λ({circumflexover (ξ)},1),

$\begin{matrix}{\mspace{79mu}{{{\hat{\xi} = \left( {{\hat{\xi}}_{1},{\hat{\xi}}_{2}} \right)},{{\hat{\xi}}_{j} = {\xi_{j}/\lambda}},{j = 1},{2\mspace{14mu}{produces}}}{{I = {\int_{0}^{\infty}{\lambda^{4}{\int_{\hat{\Omega}}{{{A\left( {\lambda\left( {\hat{\xi},1} \right)} \right)}\left\lbrack {{{w\left( {t{{x - {y(s)}}}} \right)}\left( {{t\left( {\hat{\xi},1} \right)} \cdot {\overset{.}{y}(s)}} \right)^{2}} + {O\left( \lambda^{- 1} \right)}} \right\rbrack} \times {\mathbb{e}}^{{{\mathbb{i}\lambda}{({{H{({\hat{\xi},1})}} - {{({\hat{\xi},1})} \cdot {z{({s,t})}}}})}})}{\mathbb{d}\hat{\xi}}{\mathbb{d}\lambda}}}}}},}}} & (14)\end{matrix}$where {circumflex over (Ω)} is a small neighborhood of the origin in R².The asymptotics of the integral with respect to {circumflex over (ξ)} inequation (14) is obtained by the stationary phase method:

$\begin{matrix}{{\frac{\kappa}{\lambda}{{A\left( {\lambda\left( {{\hat{\xi}\left( {s,t} \right)},1} \right)} \right)}\left\lbrack {{{w\left( {t{{x - {y(s)}}}} \right)}\left\{ {{t\left( {{\hat{\xi}\left( {s,t} \right)},1} \right)} \cdot {\overset{.}{y}(s)}} \right\}^{2}} + {O\left( \lambda^{- 1} \right)}} \right\rbrack} \times {\mathbb{e}}^{{\mathbb{i}\lambda}{({{x_{3}^{(0)}{({s,t})}} - {z_{3}{({s,t})}}})}}},{\lambda->\infty},{where}} & (15) \\{\mspace{79mu}{{{x_{3}^{(0)}\left( {s,t} \right)}\text{:}} = {H_{3}\left( {{\hat{\xi}\left( {s,t} \right)},1} \right)}}} & (16)\end{matrix}$and {circumflex over (ξ)}(s,t) is obtained by solvingH _(j)({circumflex over (ξ)},1)=z _(j)(s,t),j=1,2.  (17)In equation (15) and everywhere below κ denotes various non-zeroconstants. In equation (15) the critical point on {circumflex over (Ω)}determined by equation (17) is non-degenerate and the homogeneity ofH(ξ):H=ξ·H′={circumflex over (ξ)}·H _({circumflex over (ξ)}) +H ₃={circumflex over (ξ)}·{circumflex over (z)}+H ₃,  (18)sox ₃ ⁽⁰⁾(s,t)=H({circumflex over (ξ)}(s,t),1)−{circumflex over(ξ)}(s,t)·{circumflex over (z)}(s,t),  (19)where {circumflex over (z)}=(z₁,z₂). Substituting equations (15) and(14) into equation (13) produces:

$\begin{matrix}{{{{g(x)} - {\kappa{\int_{0}^{\infty}{\lambda^{3}{\int_{R^{2}}{{\varphi_{1}\left( {s,x} \right)}{A\left( {\lambda\left( {{\hat{\xi}\left( {s,t} \right)},1} \right)} \right)} \times \left\lbrack {{{w\left( {t{{x - {y(s)}}}} \right)}\left\{ {{t\left( {{\hat{\xi}\left( {s,t} \right)},1} \right)} \cdot {\overset{.}{y}(s)}} \right\}^{2}} + {O\left( \lambda^{- 1} \right)}} \right\rbrack \times {\mathbb{e}}^{{\mathbb{i}\lambda}\;{a{({s,t})}}}{\mathbb{d}t}{\mathbb{d}s}{\mathbb{d}\lambda}}}}}}} \in {C_{0}^{\infty}(V)}},\mspace{85mu}{{{a\left( {s,t} \right)}\text{:}} = {{x_{3}^{(0)}\left( {s,t} \right)} - {\left( {{y_{3}(s)} + {t\left( {x_{3} - {y_{3}(s)}} \right)}} \right).}}}} & (20)\end{matrix}$From equations (16) and (17),x ₃ ⁽⁰⁾(s,t)=h(ŷ(s)+t({circumflex over (x)}−ŷ(s))),  (21)where z₃=h(z₁,z₂) is the local equation of S in a neighborhood of z⁽⁰⁾.From equation (12),h ₁({circumflex over (z)} ⁽⁰⁾)=h ₂({circumflex over (z)} ⁽⁰⁾)=0.  (22)Assume first x=x⁽⁰⁾ in (20). From equation (20) minus equation (22), α=0andα_(s)=(h ₁ {dot over (y)} ₁ +h ₂ {dot over (y)} ₂ −{dot over (y)}₃)(1−t)=0,α_(t) =h ₁(x ₁ ⁽⁰⁾ −y ₁)+h ₂(x ₂ ⁽⁰⁾ −y ₂)+y ₃=0  (23)when s=s₀,t=t₀:=−c/(b−c). By assumption b≠c, so t₀ is defined. At thecritical point({circumflex over (ξ)}(s ₀ ,t ₀),1)·{dot over (y)}(s ₀)=0.  (24)From equation (12),

$\begin{matrix}{{{a_{ss}\left( {s_{0},t_{0}} \right)} = {{\left( \frac{b}{b - c} \right)^{2}\left( {{h_{11}{\overset{.}{y}}_{1}^{2}} + {2h_{12}{\overset{.}{y}}_{1}{\overset{.}{y}}_{2}} + {h_{22}{\overset{.}{y}}_{2}^{2}}} \right)} - {{\overset{¨}{y}}_{3}\frac{b}{b - c}}}},{{a_{st}\left( {s_{0},t_{0}} \right)} = {b\left( {{h_{11}{\overset{.}{y}}_{1}} + {h_{12}{\overset{.}{y}}_{2}}} \right)}},{{a_{tt}\left( {s_{0},t_{0}} \right)} = {h_{11}\left( {b - c} \right)}^{2}},{so}} & (25) \\{{\det\left( {a^{''}\left( {s_{0},t_{0}} \right)} \right)} = {{b\left\lbrack {{b{{\overset{.}{y}}_{2}^{2}\left( {{h_{11}h_{22}} - {2\left( h_{12} \right)^{2}}} \right)}} - {{\overset{¨}{y}}_{3}{h_{11}\left( {b - c} \right)}}} \right\rbrack}.}} & (26)\end{matrix}$Hence, in addition to b=0 (i.e., x⁽⁰⁾=z⁽⁰⁾), there is at most one morepoint on the ray L₊(s₀,z⁽⁰⁾), where det(α″) can be zero.Considering now x=(b,0,ρ) in (20) and arguing similarly to A. Katsevich,“Cone Beam Local Tomography”, SIAM Journal on Applied Mathematics(1999), pp. 2224-2246, gives:

$\begin{matrix}{{{s_{0}(\rho)} = {s_{0} + {O(\rho)}}},{{t_{0}(\rho)} = {t_{0} + {O(\rho)}}},{{a\left( {{s_{0}(\rho)},{t_{0}(\rho)}} \right)} = {{\frac{c}{b - c}\rho} + {O(\rho)}}},} & (27)\end{matrix}$where all O(ρ) are smooth functions of b and ρ. Consequently,({circumflex over (ξ)}(s ₀(ρ),t ₀(ρ)),1)·{dot over (y)}(s₀(ρ))=O(ρ).  (28)Using equations (27) and (28), we get under the assumption det(α″)≠0:

$\begin{matrix}{{{g\left( {b,0,\rho} \right)} = {\int_{0}^{\infty}{{\Phi\left( {{\lambda;b},\rho} \right)}{\mathbb{e}}^{{\mathbb{i}\lambda}\frac{c}{b - c}{\rho{({1 + {O{(\rho)}}})}}}{\mathbb{d}\lambda}}}},{{\Phi\left( {{\lambda;b},\rho} \right)} = {O\left( \lambda^{m + 1} \right)}},{\lambda->\infty},} & (29)\end{matrix}$where Φ admits an asymptotic expansion in decreasing powers of λ as λ→∞.The coefficients of the expansion are C^(∞) functions of b and ρ and theexpansion can be differentiated with respect to b, ρ and λ. Similarly,O(ρ) in equation (29) is a smooth function of b and ρ.

The result analogous to equation (29) for the LT function proposed in A.K. Louis and P. Maass, “Contour reconstruction in 3-D X-ray CT”, IEEETransactions on Medical Imaging, vol. 12 (1993), pp. 764-769 (see A.Katsevich, “Cone beam local tomography”, SIAM J. on Applied Mathematics,vol. 59 (1999) pp. 2224-2246) is

$\begin{matrix}{{{g_{\Lambda}\left( {b,0,\rho} \right)} = {\int_{0}^{\infty}{{\Phi\left( {{\lambda;b},\rho} \right)}{\mathbb{e}}^{{\mathbb{i}\lambda}\frac{c}{b - c}{\rho{({1 + {O{(\rho)}}})}}}{\mathbb{d}\lambda}}}},{{\Phi\left( {{\lambda;b},\rho} \right)} = {O\left( \lambda^{m + 2} \right)}},{\lambda->{\infty.}}} & (30)\end{matrix}$Equation (30) follows immediately from equations (6.10), (6.11) in A.Katsevich, “Cone beam local tomography”, SIAM J. on Applied Mathematics,vol. 59 (1999) pp. 2224-2246, if one recalls that the Fourier transformof a function in R³, which is discontinuous across a smooth surface Swith nonzero Gaussian curvature, is of order O(|ξ|⁻²) as |ξ|→∞.

Comparing equations (29) and (30) it is shown that the function inequation (2) has non-local artifacts, which are an order of magnitudesmaller. Because of the direction along which the derivative is computedin equation (2), the leading term in the expansion of Φ of orderO(λ^(m+2)) disappears and the artifact is reduced.

Dynamic Local Tomography:

Consider now the case when ƒ changes with time. The parameter s in thedefinition of the curve C is regarded as time, it is assumed that ƒ isof the form ƒ(s,x), and the cone beam data are

$\begin{matrix}{{{D_{f}\left( {s,\alpha} \right)} = {\int_{0}^{\infty}{{f\left( {s,{{y(s)} + {\alpha\; t}}} \right)}{\mathbb{d}t}}}},{s \in I},{\alpha \in {S^{2}.}}} & (31)\end{matrix}$In the dynamic case the LT function will still be defined by equation(2). The only difference is that the data are now given by equation (31)instead of equation (1).

Even though in practice ƒ(s,z) does not vanish when s∉I, to simplify thenotation assume ƒ∈C₀ ^(∞)(I×V). Clearly, this does not result in anyloss of generality, because an interval I′⊃I sufficiently close to I isselected, C is extended appropriately, and ƒ is multiplied by a cut-offχ∈C₀ ^(∞)(I′), such that χ=1 on I. Since φ∈C₀ ^(∞)(I×V), replacing ƒ byχƒ will not affect g. The same argument applies when ƒ is adistribution.

First, the operator B defined by equations (31) and (2) extends to a mapE′(I×V)→E′(V).

$\begin{matrix}{{{{For}\mspace{14mu}\psi} \in {{C^{\infty}(V)}\left( {{Bf},\psi} \right)}} = {{\int_{V}{{\psi(x)}{\int_{I}{{\varphi\left( {s,x} \right)} \times {\frac{\partial^{2}}{\partial q^{2}}\left\lbrack {{{x - {y(q)}}}{\int_{0}^{\infty}{{w\left( {t{{x - {y(q)}}}} \right)}{f\left( {s,{{y(s)} + {t\left( {x - {y(q)}} \right)}}} \right)}{\mathbb{d}t}}}} \right\rbrack}}}}}❘_{q = s}{{\mathbb{d}s}{{\mathbb{d}x}.}}}} & (32)\end{matrix}$Changing variables x→z=y(s)+t(x−y(q)), after simple transformations

$\begin{matrix}{\left( {{Bf},\psi} \right) = {\int_{I}{\int_{V}{\left\lbrack {{{{z - {y(s)}}}{w\left( {t{{z - {y(s)}}}} \right)}\frac{\partial^{2}}{\partial q^{2}}\left( {\int_{0}^{\infty}{{\varphi\left( {s,{x( \cdot )}} \right)}{\psi\left( {s,{x( \cdot )}} \right)}t^{- 4}{\mathbb{d}t}}} \right)}❘_{q = s}} \right\rbrack \times {f\left( {s,z} \right)}{\mathbb{d}z}{{\mathbb{d}s}.{where}}}}}} & (33) \\{{{x( \cdot )}\text{:} = {x\left( {z,s,q,t} \right)}} = {{t^{- 1}\left( {z - {y(s)}} \right)} + {{y(q)}.}}} & (34)\end{matrix}$

Since φ∈C₀ ^(∞)(I×V), t is bounded away from zero in equation (33), theexpression in brackets is a C^(∞)(I×V) function, and the desiredassertion follows immediately by duality.

Next the relation between the wave fronts of ƒ and Bƒ is investigatedwithout considering the most general situation, but concentrating on thepractically relevant case when ƒ is a conormal distribution, whichdepends smoothly on s. An arbitrary (s₀,z⁽⁰⁾,η⁽⁰⁾)∈I×(T*V\0) isselected. Because of linearity, the following assumptions are made aboutƒ without loosing any generality.

-   -   1. supp ƒ belongs to a small neighborhood of (s₀,z⁽⁰⁾);    -   2. For a sufficiently small open cone Ω        η⁽⁰⁾ the Fourier transform of ƒ(s,z) with respect to z        satisfies:        {tilde over (ƒ)}(s,η)=A(s,η)e ^(iH(s,η)) +A ₁(s,η).  (35)    -   3. Here        -   A(s,η) is smooth, identically vanishes on R³, Ω, and            (∂/∂s)^(k)A(s,η)∈S^(m)(R³) uniformly in s for all k≧0;        -   H(s,η) is smooth (away from ξ=0), real-valued, and            homogeneous of degree one in η;        -   A₁(s,η)∈C^(∞) and (∂/∂s)^(k)A₁(s,η)=o(|η|^(−N)),|η|→∞,            uniformly in s for all k≧0 and N≧0.

When ƒ is as described above and WF(g) is a subset of all (x,ξ)∈T*V\0which satisfyy(s)+t(x−y(s))=H _(η)(s,η)η·{dot over (y)}(s)(1−t)=H _(s)(s,η)η·(x−y(s))=0ξ=tη  (36)for some (s,x)∈supp φ and η∈Ω.

Clearly the derivative ∂²/∂q² in equation (2) can be omitted whenfinding WF(g). From equation (35) the integral is studied

$\begin{matrix}{\int_{R^{3}}{\int_{\Omega}{\int_{I}{\int_{0}^{\infty}{{\varphi_{1}\left( {s,t,x} \right)}{A\left( {s,\eta} \right)}{\mathbb{e}}^{{\mathbb{i}}{\lbrack{{H{({s,\eta})}} - {\eta \cdot {({{y{(s)}} + {t{({x - {y{(s)}}})}}})}} + {\xi \cdot x}}\rbrack}}{\mathbb{d}t}{\mathbb{d}s}{\mathbb{d}\eta}{\mathbb{d}x}}}}}} & (37)\end{matrix}$for some smooth φ₁ with compact support in the domain t>0.Differentiating the expression in brackets in equation (37) produces[·]_(η) =H _(η)(s,η)−(y(s)+t(x−y(s)))[·]_(s) =H _(s)(s,η)−η·{dot over (y)}(s)(1−t)[·]_(t)=−η·(x−y(s))[·]_(x) =ξ−tη.  (38)

In an embodiment, a shift between the location of a reconstruction pointshown on the reconstructed image and the actual location of the point ofthe useful feature. The step of determining the error of the movingobject is especially useful when reconstructing an image of a movingobject to correct for the location of the useful features of the actualimage. Suppose first that η·{dot over (y)}(s)≠0 in equation (36). Thento (z,η)∈WF(ƒ(s,·)) such that y(s)∈Π(z,η)∩C there corresponds(x,ξ)∈WF(g), where

$\begin{matrix}{{x = {z + {\left( {z - {y(s)}} \right)\left( {t^{- 1} - 1} \right)}}},{\xi = {t\;\eta}},{t = {1 - {\frac{H_{s}\left( {s,\eta} \right)}{\eta \cdot {\overset{.}{y}(s)}}.}}}} & (39)\end{matrix}$If η·{dot over (y)}(s)=0 and H_(s)(s,η)=0, then to (z,η)∈WF(ƒ(s,·)) suchthat y(s)∈Π(z,η)∩C there corresponds (x,ξ)∈WF(g), wherex=y(s)+(z−y(s))/t,ξ=tη,t>0.  (40)

Thus, similarly to the static case, WF(g) consists of usefulsingularities described by equation (39), and of non-local artifacts,which are described by equation (40). However, as opposed to the staticcase, the useful singularities of g in general no longer coincide withWF(ƒ). For example, the size of the shift between singsupp g andsingsupp ƒ depends on the quantity H_(s)(s,η)/(η·{dot over (y)}(s)).Knowing this quantity can lead to the increased accuracy of locatingsingsupp ƒ from singsupp g. If the motion of the object is approximatelyknown, one can locate the singularities of g (represented as(x,ξ)∈WF(g)), then estimate the quantity H_(s)(s,η)/(η·{dot over(y)}(s)) based upon a priory information and equation (36), and then useequation (39) to find singsupp ƒ (represented by z in that equation).

The conditions leading to the appearance of the non-local artifact aresomewhat different compared with the static case as well. The rayL₊(s,z)⊂g only when the additional condition H_(s)(s,η)=0 holds. Notethat when H_(s)(s,η)≡0, equation (36) transforms to equation (5).

Next consider the quantitative relationship between the singularities ofƒ and g. First, lemma 3 is needed.

Suppose Π(z⁽⁰⁾,η⁽⁰⁾)∩C=y(s₀) and η⁽⁰⁾·{dot over (y)}(s₀)≠0, wherez⁽⁰⁾=H_(η)(s₀,η⁽⁰⁾). Find (x⁽⁰⁾,ξ⁽⁰⁾) by solving equation (36). Let Ω bea sufficiently small open cone containing ξ⁽⁰⁾. Solving the systemequation (36) for s,t,x, and η in terms of ξ∈Ω determines C^(∞)(Ω)functions s(ξ),t(ξ),x(ξ), and η(ξ). s,t,x are homogeneous of order 0,and η is homogeneous of order 1.

Choose the coordinate system in which the x₃-coordinate is along η⁽⁰⁾.The homogeneity of H implies ∂²/(∂η_(j)∂η₃)H(s,η)|_(η=η) ₍₀₎ =0,j=1,2,3. Using equation (38), the matrix of the second derivativesevaluated at (s₀,x⁽⁰⁾,η⁽⁰⁾,ξ⁽⁰⁾,t), where t is determined from thesecond equation in equation (36), becomes

$\begin{matrix}{{D\text{:}} = {{\begin{matrix}H_{11}^{''} & H_{12}^{''} & 0 & {H_{s\; 1}^{''} - {{\overset{.}{y}}_{1}\left( {1 - t} \right)}} & {- \Delta_{1}} & {- t} & 0 & 0 \\H_{21}^{''} & H_{22}^{''} & 0 & {H_{s\; 2}^{''} - {{\overset{.}{y}}_{2}\left( {1 - t} \right)}} & {- \Delta_{2}} & 0 & {- t} & 0 \\0 & 0 & 0 & {H_{s\; 3}^{''} - {{\overset{.}{y}}_{3}\left( {1 - t} \right)}} & 0 & 0 & 0 & {- t} \\\; & \left. {{\overset{.}{y}}_{2}\left( {1 - t} \right)} \right\rbrack & \left. {{\overset{.}{y}}_{3}\left( {1 - t} \right)} \right\rbrack & \left. {\eta_{3}{{\overset{¨}{y}}_{3}\left( {1 - t} \right)}} \right\rbrack & {\eta_{3}{\overset{.}{y}}_{3}} & 0 & 0 & 0 \\{- \Delta_{1}} & {- \Delta_{2}} & 0 & {\eta_{3}{\overset{.}{y}}_{3}} & 0 & 0 & 0 & {- \eta_{3}} \\{- t} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & {- t} & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & {- t} & 0 & {- \eta_{3}} & 0 & 0 & 0\end{matrix}}.}} & (41)\end{matrix}$Here Δ_(j)=x_(j)−y_(j)(s), H″_(jk)=∂²/(∂η_(j)∂η_(k))H(s,η),H″_(sj)=∂²/(∂s∂η_(j))H(s,η), j,k=1,2,3. In equation (41), Δ₃=0. Aftersimple transformations,D=t ⁴η₃ ²({dot over (y)} ₃ −H _(s3))².  (42)Because of the homogeneity of H, H_(s)=ξ·H_(sξ). Using equations (36)and (42) producesD=t ⁴(η·{dot over (y)}−η·H _(sη))² =t ⁴(η·{dot over (y)}−H _(s))² =t⁶(η·{dot over (y)})²>0.  (43)This proves that the functions s(ξ),t(ξ),x(ξ), and η(ξ) are C^(∞) nearξ⁽⁰⁾. The statement about homogeneity follows from equation (36) and thehomogeneity of H. Therefore, the functions are C^(∞) inside the cone Ω.

Using G(ξ):=ξ·x(ξ), the Lagrangian submanifold Λ of a conic neighborhoodof (x⁽⁰⁾,ξ⁽⁰⁾) is defined as: Λ={(G′(ξ),ξ):ξ∈Ω}.

Again using ƒ as previously described, suppose Π(z⁽⁰⁾,η⁽⁰⁾)∩C=y(s₀) andη⁽⁰⁾·{dot over (y)}(s₀)≠0. (x⁽⁰⁾,ξ⁽⁰⁾) is found by solving equation(36). Then equation (2) defines a distribution g∈I^(m+7/4)(V,Λ), and for|ξ|>1:

$\begin{matrix}{{{{{\mathbb{e}}^{{\mathbb{i}}\;{G{(\xi)}}}{\overset{\_}{g}(\xi)}} + {2\pi\frac{\varphi\left( {s,x} \right){{x - {y(s)}}}}{t}{{\eta \cdot {\overset{.}{y}(s)}}}{A\left( {s,\eta} \right)}}} \in {S^{m}\left( R^{3} \right)}},} & (44)\end{matrix}$where s,t,x, and η are the functions of ξ obtained by solving the systemequation (36), and A(s,η) is interpreted as 0 if no solution exists.e^(−iG(ξ)){tilde over (g)}(ξ) is polyhomogeneous if A is.

Similarly to equations (8) and (9), it is found that

$\begin{matrix}{{\overset{\sim}{g}({\sigma\xi})} = {{- \frac{1}{\left( {2\pi} \right)^{3}}}{\int_{V}{\int_{R^{3}}{\int_{I}{\int_{0}^{\infty}{{{\varphi\left( {s,x} \right)}\left\lbrack {{{w\left( {t{{x - {y(s)}}}} \right)}{{x - {y(s)}}}\left( {t\;{\eta \cdot {y(s)}}} \right)^{2}} + {O\left( {\eta } \right)}} \right\rbrack} \times {A\left( {s,\eta} \right)}{\mathbb{e}}^{{\mathbb{i}}{\lbrack{{H{({s,\eta})}} - {\eta \cdot {({{y{(s)}} + {t{({x - {y{(q)}}})}}})}} + {{\sigma\xi} \cdot x}}\rbrack}}{\mathbb{d}t}{\mathbb{d}s}{\mathbb{d}\eta}{{\mathbb{d}x}.}}}}}}}} & (45)\end{matrix}$Changing variables η→ση and factoring out σ, the same expression as theone in brackets in equation (37). Equation (43) implies that thesignature of the matrix in equation (41) is a multiple of 4. Applyingthe same technique as in described in “The analysis of linear partialdifferential operators”, to the integral with respect to η in equation(45) and then using the stationary phase method.

This method allows us to generalize the notion of visible singularitiesto the dynamic case. The singularity (z,η)∈WF(ƒ(s,·)) is visible, i.e.can be stably detected from the data, if (z−y(s))·η=0 and η·{dot over(y)}(s)≠0. Comparing the results with equations (8) and (11), it isconcluded that in both the static and dynamic cases the visiblesingularities of ƒ are reconstructed with the same resolution. Inparticular, there is no smearing due to motion in the dynamic case.

As previously described, non-local artifacts arise in the dynamic caseand the strength of these artifacts is the same as in the static case.

(z⁽⁰⁾=H_(η)(s₀,η⁽⁰⁾),η⁽⁰⁾)∈WF(ƒ) is selected such that the planeΠ(z⁽⁰⁾,η⁽⁰⁾) is tangent to C at y(s₀) and H_(s)(s₀,η⁽⁰⁾)=0. Fixx⁽⁰⁾∈L₊(s₀,z⁽⁰⁾), x⁽⁰⁾∉{z⁽⁰⁾,y(s₀)} and (s₀,x⁽⁰⁾)∈suppφ. It is againassumed that ƒ is of the form previously described and the Gaussiancurvature of the surface z(η)=H′_(η)(s₀,η),η∈Ω, is nonzero near z⁽⁰⁾.Subtracting equation (21) from equation (12) produces the sameexpression as in equation (20), but the function x₃ ⁽⁰⁾ is now given byx ₃ ⁽⁰⁾(s,t)=h(s,ŷ(s)+t({circumflex over (x)}−ŷ(s))).  (46)Here z₃=h(s,z₁,z₂) is determined by solvingH _(j)(s,({circumflex over (ξ)},1))=z _(j)(s,t),j=1,2,  (47)for {circumflex over (ξ)}={circumflex over (ξ)}(s,t) (cf. (17)), andthen substituting {circumflex over (ξ)}(s,t) into H′₃ (cf. (16) and(21)):h(s,z ₁ ,z ₂):=H ₃(s,({circumflex over (ξ)}(s,t),1)).  (48)By assumption,h ₁(s ₀ ,{circumflex over (z)} ⁽⁰⁾)=h ₂(s ₀ ,{circumflex over (z)}⁽⁰⁾)=0.  (49)Again assuming x=x⁽⁰⁾, we get α=0 andα_(s) =h _(s)+(h ₁ {dot over (y)} ₁ +h ₂ {dot over (y)} ₂ −{dot over(y)} ₃)(1−t)=0,α_(t) =h ₁(x ₁ ⁽⁰⁾ −y ₁)+h ₂(x ₂ ⁽⁰⁾ −y ₂)+y ₃=0  (50)when s=s₀,t=t₀:=−c/(b−c). The additional condition H_(s)(s₀,η⁽⁰⁾)=0 isused which implies h_(s)(s₀,{circumflex over (z)}⁽⁰⁾)=0. The analogue ofequation (25) becomes

$\begin{matrix}{{{a_{ss}\left( {s_{0},t_{0}} \right)} = {h_{ss} + {2\left( {{h_{s\; 1}{\overset{.}{y}}_{1}} + {h_{s\; 2}{\overset{.}{y}}_{2}}} \right)\frac{b}{b - c}} + {\left( \frac{b}{b - c} \right)^{2}\left( {{h_{11}{\overset{.}{y}}_{1}^{2}} + {2h_{12}{\overset{.}{y}}_{1}{\overset{.}{y}}_{2}} + {h_{22}{\overset{.}{y}}_{2}^{2}}} \right)} - {{\overset{¨}{y}}_{3}\frac{b}{b - c}}}},{{a_{st}\left( {s_{0},t_{0}} \right)} = {{h_{s\; 1}\left( {b - c} \right)} + {b\left( {{h_{11}{\overset{.}{y}}_{1}} + {h_{12}{\overset{.}{y}}_{2}}} \right)}}},{{a_{tt}\left( {s_{0},t_{0}} \right)} = {h_{11}\left( {b - c} \right)}^{2}},{so},} & (51) \\{{\det\left( {a^{''}\left( {s_{0},t_{0}} \right)} \right)} = {{\left( {b - c} \right)^{2}\left( {{h_{ss}h_{11}} - \left( h_{s\; 1} \right)^{2}} \right)} + {2{b\left( {b - c} \right)}{{\overset{.}{y}}_{2}\left( {{h_{11}h_{s\; 2}} - {h_{s\; 1}h_{12}}} \right)}} + {{b\left\lbrack {{b{{\overset{.}{y}}_{2}^{2}\left( {{h_{11}h_{22}} - {2\left( h_{12} \right)^{2}}} \right)}} - {{\overset{¨}{y}}_{3}{h_{11}\left( {b - c} \right)}}} \right\rbrack}.}}} & (52)\end{matrix}$There are at most two more points on the ray L₊(s₀,z⁽⁰⁾) (i.e., twovalues of b) where det(α″) can be zero. Thus, under the assumption thatdet(α″)≠0, the same expression as in equation (29) is derived for thenon-local artifact in the dynamic case.Numerical Implementation and Experiments:

Since the LT function in equation (2) stays the same regardless ofwhether ƒ depends on time or not, implementation of equation (2) isidentical in both the static and dynamic cases. The followingdescription describes how equation (2) is computed when the sourcetrajectory is a helix:

$\begin{matrix}{{{y_{1}(s)} = {R\;{\cos\left( {s - s_{0}} \right)}}},{{y_{2}(s)} = {R\;{\sin\left( {s - s_{0}} \right)}}},{{y_{3}(s)} = {y_{3}^{(0)} + {\frac{h}{2\pi}{s.}}}}} & (53)\end{matrix}$Here

-   -   R is the radius of the helix;    -   h is the pitch of the helix;    -   s₀ and y₃ ⁽⁰⁾ are the parameters that determine the initial        point on the helix; and    -   (y₁(s),y₂(s),y₃(s)) are the coordinates of a point on the helix.

Following the curved detector geometry described in F. Noo, J. Pack andD. Heuscher, “Exact helical reconstruction using native cone-beamGeometries”, Physics in Medicine and Biology, vol. 48, (2003), pp.3787-3818, let α and w be the conventional detector coordinates: α isthe polar angle (with the origin at the source), and w is the verticalcoordinate. Hence the cone beam data are given by D_(ƒ)(s,(α,w)).Assuming that the detector contains the axis of rotation, let DP(s)denote the surface of the detector, which corresponds to the sourcelocated at y(s). In the numerical experiments it is tacitly assumed thatthe detector is large enough, e.g. that it contains at least the Tamwindow. In general, however, this is not required for LT. Let (α(q,s,x),w(q,s,x)) be the coordinates of the point where the rayz=y(s)+tβ(q,x),t>0, intersects DP(s). Similarly, let (α(s,x), w(s,x)) bethe coordinates of the projection of x onto DP(s). Clearly,α(q,s,x)|_(q=s)=α(s,x),w(q,s,x)|_(q=s) =w(s,x).  (54)Using the chain rule:

$\begin{matrix}\begin{matrix}{{g(x)} = {{\int_{I}{{\varphi\left( {s,x} \right)}\frac{\partial^{2}}{\partial q^{2}}{D_{f}\left( {s,\left( {{\alpha\left( {q,s,x} \right)},{w\left( {q,s,x} \right)}} \right)} \right)}}}❘_{q = s}{\mathbb{d}s}}} \\{{{\cong {\int_{I}{{\varphi\left( {s,x} \right)}\left( {{c_{\alpha}^{2}\frac{\partial^{2}}{\partial\alpha^{2}}} + {2c_{\alpha}c_{w}\frac{\partial^{2}}{{\partial\alpha}{\partial w}}} + {c_{w}^{2}\frac{\partial^{2}}{\partial w^{2}}}} \right){D_{f}\left( {s,\left( {\alpha,w} \right)} \right)}}}}❘_{\underset{w = {w{({s,x})}}}{a = {a{({s,x})}}}}{\mathbb{d}s}},}\end{matrix} & (55)\end{matrix}$is obtained from equation (2) where

$\begin{matrix}{{c_{\alpha} = {\frac{\partial{\alpha\left( {q,s,x} \right)}}{\partial q}❘_{q = s}}},\mspace{79mu}{c_{w} = {\frac{\partial{w\left( {q,s,x} \right)}}{\partial q}❘_{q = s}}},} & (56)\end{matrix}$and ≅ denotes equality up to the leading singular term. The termsomitted on the right in equation (55) contain, at most, the first orderderivatives of D_(ƒ), so that the results previously obtained apply tothe new function.

Using the approach described in id., for α(q,s,x) and w(q,s,x):

$\begin{matrix}{{{{\tan\;{\alpha\left( {q,s,x} \right)}} = \frac{{{y_{1}(s)}\left( {x_{2} - {y_{2}(q)}} \right)} - {{y_{2}(s)}\left( {x_{1} - {y_{1}(q)}} \right)}}{{{y_{1}(s)}\left( {x_{1} - {y_{1}(q)}} \right)} - {{y_{2}(s)}\left( {x_{2} - {y_{2}(q)}} \right)}}},{{w\left( {q,s,x} \right)} = {R\;{\frac{x_{3} - {y_{3}(q)}}{\sqrt{\left( {x_{1} - {y_{1}(q)}} \right)^{2} + \left( {x_{2} - {y_{2}(q)}} \right)^{2}}}.{Recall}}\mspace{14mu}{that}\text{:}}}}\mspace{574mu}} & (57) \\{{{V\left( {s,x} \right)}:={R - \frac{{x_{1}{y_{1}(s)}} + {x_{2}{y_{2}(s)}}}{R}}},{{\tan\;{\alpha\left( {s,x} \right)}} = \frac{{x_{2}{y_{1}(s)}} + {x_{1}{y_{2}(s)}}}{{RV}\left( {s,x} \right)}},{{w\left( {s,x} \right)} = {R\;\cos\;{\alpha\left( {s,x} \right)}{\frac{x_{3} - {y_{3}(s)}}{V\left( {s,x} \right)}.}}}} & (58)\end{matrix}$Differentiating equation (57) with respect to q and using equations (58)and (54), after some transformations, it is found that

$\begin{matrix}{{c_{\alpha} = {- \frac{R\;\cos^{2}{\alpha\left( {s,x} \right)}}{V\left( {s,x} \right)}}},{c_{w} = {\frac{R\;\cos\;{\alpha\left( {s,x} \right)}}{V\left( {s,x} \right)}{\left( {{{w\left( {s,x} \right)}\sin\;{\alpha\left( {s,x} \right)}} - \frac{h}{2\pi}} \right).\square^{*}}}}} & (59)\end{matrix}$

In general, there is a significant flexibility in selecting the functionφ(s,x). By choosing the appropriate φ(s,x), a LT can be adopted forpractically any scan geometry, including truncated projections,incomplete or segmented source trajectory, etc. A possible candidate,which guarantees the reconstruction of all visible singularities and isefficient from the computational point of view is determined. Letw=w_(bot)(α) and w=w_(top)(α) be the equations of the top and bottomboundaries of the Tam window on the detector. Choose any χ∈C₀ ^(∞)(R)such that χ≡1 on [0,1], and define

$\begin{matrix}{{\varphi\left( {s,x} \right)} = {{\chi\left( \frac{{w_{top}(\alpha)} - w}{{w_{top}(\alpha)} - {w_{bot}(\alpha)}} \right)}❘_{\underset{w = {w{({s,x})}}}{a = {a{({s,x})}}}}.}} & (60)\end{matrix}$By construction, φ(s,x)>0 whenever s∈I_(π)(x) (the π-interval of x).This property guarantees that the visible singularities of ƒ aredetected (see equations (11) and (44)). The advantage of equation (60)is that φ(s,x) only depends on the projection of x onto DP(s), so thefunction is written in the form φ(α(s,x), w(s,x)). Equation (59) is usedto rewrite equation (55) in the form

$\begin{matrix}{\mspace{79mu}{{{g(x)} = {\int_{I}{\frac{R^{2}}{V^{2}\left( {s,x} \right)}{\Phi\left( {{s;{\alpha\left( {s,x} \right)}},{w\left( {s,x} \right)}} \right)}{\mathbb{d}s}}}},{{\Phi\left( {{s;\alpha},w} \right)} = {{\varphi\left( {\alpha,w} \right)}\cos^{2}{\alpha\left\lbrack {{\cos^{2}\alpha\frac{\partial^{2}}{\partial\alpha^{2}}} - {2\cos\;{\alpha\left( {{w\;\sin\;\alpha} - {h/\left( {2\pi} \right)}} \right)}\frac{\partial^{2}}{{\partial\alpha}{\partial w}}} + {\left( {{w\;\sin\;\alpha} - {h/\left( {2\pi} \right)}} \right)^{2}\frac{\partial^{2}}{\partial w^{2}}}} \right\rbrack}{{D_{f}\left( {s,\left( {\alpha,w} \right)} \right)}.}}}}} & (61)\end{matrix}$As is seen, equation (61) admits a very efficientfiltered-backprojection implementation. Moreover, the filtering stepconsists only of computing derivatives, i.e. is local.

For comparison purposes the LT function is computed which is analogousto the one described in A. K. Louis and P. Maass, “Contourreconstruction in 3-D X-ray CT”, IEEE Transaction on Medical Imaging,vol. 12 (1993) pp. 764-769 as:

$\begin{matrix}\begin{matrix}{{g_{\Lambda}(x)} = {\int_{I}{{\varphi\left( {s,x} \right)}\Delta_{x}{D_{f}\left( {s,\left( {{\alpha\left( {s,x} \right)},{w\left( {s,x} \right)}} \right)} \right)}{\mathbb{d}s}}}} \\{\cong {\int_{I}{{\varphi\left( {s,x} \right)}\left( {{{{\nabla_{x}{\alpha\left( {s,x} \right)}}}^{2}\frac{\partial^{2}}{\partial\alpha^{2}}} +} \right.}}} \\{{2\left( {{\nabla_{x}{\alpha\left( {s,x} \right)}} \cdot {\nabla_{x}{w\left( {s,x} \right)}}} \right)\frac{\partial^{2}}{{\partial\alpha}{\partial w}}} +} \\{{{\left. {{{\nabla_{x}{w\left( {s,x} \right)}}}^{2}\frac{\partial^{2}}{\partial w^{2}}} \right){D_{f}\left( {s,\left( {\alpha,w} \right)} \right)}}❘_{\underset{w = {w{({s,x})}}}{a = {a{({s,x})}}}}{\mathbb{d}s}},}\end{matrix} & (62)\end{matrix}$Using equation (58) produces

$\begin{matrix}{{{\nabla_{x}{\alpha\left( {s,x} \right)}} = {{\frac{\cos\;\alpha}{RV}\left( {{{y_{1}\sin\;\alpha} - {y_{2}\cos\;\alpha}},{{y_{1}\cos\;\alpha} + {y_{2}\sin\;\alpha}},0} \right)}❘_{\underset{V = {V{({s,x})}}}{a = {a{({s,x})}}}}}},{{\nabla_{x}{w\left( {s,x} \right)}} = {{\frac{\cos\;\alpha}{RV}\text{(}w\text{(}y_{1}\cos\;\alpha} + {y_{2}\sin\;\alpha\text{)}}}},{{w\text{(}} - {y_{1}\sin\;\alpha} + {y_{2}\cos\;\alpha\text{)}}},{{R^{2}\text{)}}❘_{\underset{\underset{w = {w{({s,x})}}}{V = {V{({s,x})}}}}{a = {a{({s,x})}}}}.}} & (63)\end{matrix}$After simple calculations, equation (63) implies

$\begin{matrix}{{{{\nabla_{x}\alpha}}^{2} = \left( \frac{\cos\;\alpha}{V} \right)^{2}},{{{\nabla_{x}\alpha} \cdot {\nabla_{x}w}} = 0},{{{\nabla_{x}w}}^{2} = {\left( \frac{\cos\;\alpha}{V} \right)^{2}{\left( {R^{2} + w^{2}} \right).}}}} & (64)\end{matrix}$Finally, substitution of equation (64) into equation (62) yields

$\begin{matrix}{{{g_{\Lambda}(x)} = {\int_{I}{\frac{1}{V^{2}\left( {s,x} \right)}{\Phi_{\Lambda}\left( {{s;{\alpha\left( {s,x} \right)}},{w\left( {s,x} \right)}} \right)}{\mathbb{d}s}}}},{{\Phi_{\Lambda}\left( {{s;\alpha},w} \right)} = {{\varphi\left( {\alpha,w} \right)}\cos^{2}{\alpha\left\lbrack {\frac{\partial^{2}}{\partial\alpha^{2}} + {\left( {R^{2} + w^{2}} \right)\frac{\partial^{2}}{\partial w^{2}}}} \right\rbrack}{{D_{f}\left( {s,\left( {\alpha,w} \right)} \right)}.}}}} & (65)\end{matrix}$

FIG. 1 is a flow diagram of the steps for the cone beam local tomographyaccording to the present invention. In step 10 the current cone beamprojection is loaded into computer memory with the assumption that theCB projection loaded into computer memory corresponds to the x-raysource located at y(s). Using the CB projection data D_(ƒ)(y(s),(α,w))in step 20, the combination of derivatives that will result insuppression of the artifacts of the image is determined, wherein thecombination of derivatives corresponds to the differentiation of thedata along the direction tangent to the curve. Based on the determinedcombination direction, the cone beam projection data is differentiated.In this example, the derivatives

${\frac{\partial^{2}}{\partial\alpha^{2}}{D_{f}\left( {s,\left( {\alpha_{i},w_{j}} \right)} \right)}},{\frac{\partial^{2}}{{\partial\alpha}{\partial w}}{D_{f}\left( {s,\left( {\alpha_{i},w_{j}} \right)} \right)}},{\frac{\partial^{2}}{\partial w^{2}}{D_{f}\left( {s,\left( {\alpha_{i},w_{j}} \right)} \right)}}$are found for all (α_(i),w_(j)) on the detector. The computation can beaccomplished using interpolation and finite differences. Thus, step 20includes finding the combination of derivatives that will result inreduced artifacts and then computing the derivatives, in other words, itis local.

In step 30, the derivative results are processed to compute thedetermined combination of the derivatives. In this example, thederivatives are multiplied by the quantities cos² α, −2 cos α(w sinα−h/(2π)), and (w sin α−h/(2π))², respectively, for all (α_(i), w_(j))on the detector. Then the results are multiplied by φ(α_(i), w_(j))cos²α_(i) to produce the function Φ(s; α_(i), w_(j)) according to equation(61), where the function φ(α, w) is determined using equation (60) andthe relation φ(α(s,x), w(s,x))=φ(s,x). By construction, φ(s, x)>0whenever s∈I_(π)(x) (the π-interval of x). This property guarantees thatthe visible singularities of ƒ are detected. The advantage of equation(60) is that φ(s, x) only depends on the projection of x onto DP(s), sothe function is written in the form φ(α(s,x), w(s,x)).

In backprojection step 40, for each reconstruction point x, the datafound in step 30 is back-projected according to the first line ofequation (61). FIG. 2 is a seven substep flow diagram forbackprojection, which corresponds to step 40 of FIG. 1. In step 41, areconstruction point x is fixed, which represents a point inside thepatient where it is required to reconstruct the image. Thereconstruction point x is projected onto the detector in step 42 bycomputing α(s,x), w(s,x) by formulas given in equation (58). In step 43it is determined if φ(α(s,x), w(s,x))=0. If φ(α(s,x), w(s,x))=0, thenthe preprocessed CB data is not used for image reconstruction at x andstep 41 is repeated for a different reconstruction point. If φ(α(s,x),w(s,x))≠0, then the filtered CB data affects the image at x and steps44-47 are performed.

In step 44, the value Φ(s;α(s,x), w(s,x)) is estimated by interpolationfrom a few values Φ(s;α_(i),w_(j)) for α_(i),w_(j) close to α(s,x),w(s,x). Φ(s;α(s,x), w(s,x)) is multiplied by

$\frac{R^{2}}{V^{2}\left( {s,x} \right)}$in step 45, where V(s,x) is computed by equation (58). The result ofstep 45 is added to the image being reconstructed at the point x in step46 according to a pre-selected scheme (for example, the Trapezoidalscheme) for approximate evaluation of the integral in equation (61). Instep 47, go to step 41 and select a different reconstruction point x.After repeating steps 41-47 for each reconstruction point inside thearea of interest, the process advances to step 50.

Referring back to FIG. 1, steps 10-40 are repeated until the scan isfinished or when the image reconstruction process is complete for therequired points of the area of interest to reconstruct the usefulfeatures. At step 50, the process returns to step 10 shown in FIG. 1 toload the next CB projection data into computer memory. The image can bedisplayed at all reconstruction points x for which the imagereconstruction process has been completed, that is, when subsequent CBprojections are not necessary for reconstructing the image at thosepoints. The CB projection that have already been processed is no longerneeded for image reconstruction and is discarded from the computermemory. The algorithm concludes when the scan is finished or the imagereconstruction process has completed at all the required points.

For LT function testing, experiments were performed. Table 1 providessimulation and reconstruction parameters used for the experiments.

TABLE 1 Simulation and reconstruction parameters R (radius of thespiral) 600 mm h (pitch of the spiral)  66 mm detector pixel size 10⁻³rad × 0.6 mm (as projected to isocenter) number of detector rows 128number of detectors per row 950 number of source positions per rotation1000  voxel size in each direction  1.0 mm

FIGS. 3 a and 3 b show the results of reconstructing the conventionalclock phantom as described in A. Katsevich, Samit Basu, and Jiang Hsieh,“Exact filtered backprojection reconstruction for dynamic pitch helicalcone beam computed tomography”, Physics in Medicine and Biology, vol. 49(2004) pp. 3089-3103. The original clock phantom, which is slightlydifferent from the one used here is described in H. Turbell and P. E.Danielsson, “Helical cone beam tomography”, Int. Journal of ImagingSystem and Technology, vol. 11 (2000), pp. 91-100.

FIG. 4 a shows the LT function g according to the present invention andFIG. 4 b shows the LT function g_(Λ) according to the prior art. FIG. 4a shows the results of reconstructing the phantom, which consists ofthree elongated ellipsoids with half-axes 100×100×10 (all the sizes andcoordinates are given in mm). The locations of their centers are(0,0,−40), (0,0,0), and (0,0,40), and the cross-section|x₁|≦150,x₂=0,|x₃|≦60 is shown.

FIGS. 5 a and 5 b show the results for another phantom according toanother experiment. The phantom in this embodiment consists of twoidentical balls having radius 40 centered at approximately (−80,0,0) andapproximately (80,0,0), respectively. As is seen from the figures, thenon-local artifacts in g are weaker than those in g_(Λ). Reconstructionof the clock phantom with the cross-section |x₁|≦255.5,|x₂|≦255.5,x₃=0is shown in FIG. 3 a. FIG. 3 a shows the LT function g according to thepresent invention and FIG. 3 b shows the LT function g_(Λ) according tothe prior art.

Reconstructions in the dynamic case are shown in FIGS. 6 and 7 whereinthe phantom is represented by a single ellipsoid

$\begin{matrix}{{{\left( \frac{x_{1} - {x_{1}(s)}}{a(s)} \right)^{2} + \left( \frac{x_{2} - {x_{2}(s)}}{b(s)} \right)^{2} + \left( \frac{x_{3} - {x_{3}(s)}}{c(s)} \right)} \leq 1},} & (66)\end{matrix}$wherein the half-axes (a(s), b(s), c(s)) and the center(x₁(s),x₂(s),x₃(s))—depend on time. For FIG. 6 a:α(s)=b(s)=c(s)=50,x ₁(s)=x ₂(s)=0,x ₃(s)=10 sin(s/3).  (67)For FIG. 6 b:α(s)=b(s)=c(s)=50,x ₁(s)=10 cos(s/3),x ₂(s)=x ₃(s)=0.  (68)For FIG. 7 a:α(s)=b(s)=50(1+0.2 cos(s/3)),c(s)=50(1+0.2 sin(s/3)),x ₁(s)=x ₂(s)=x₃(s)=0.  (69)For FIG. 7 b:α(s)=b(s)=50(1+0.2 cos(s/3)),c(s)=50(1+0.2 sin(s/3)),x ₁(s)=10cos(s/3),x ₂(s)=x ₃(s)=0.  (70)The results of dynamic reconstructions confirm that the visiblesingularities are reconstructed with about the same resolution as in thestatic case, and that no additional artifacts arise because of thechanges in the phantom. As soon as some singularity of ƒ becomes“invisible” due to motion, the corresponding singularity in g ends atthat point, but this does not cause any streaks across the image.Equivalently we can say that LT is quite stable with respect toinconsistencies in the data caused by motion.

The images shown in FIGS. 6 and 7 were obtained by ignoring the dynamicnature of the phantom during reconstruction. Clearly, it istheoretically impossible to stably reconstruct an accurate 4D image of ƒ(3D+time) from the cone beam data with only three degrees of freedom. Onthe other hand, in many cases one can use additional information forimproving image quality. Consider, for example, cardiac imaging based ona circular source trajectory, which is of significant importance inmedical CT.

Using the electro cardiogram (ECG) of the patient, which is measuredconcurrently with cone beam projections, the data corresponding to afixed cardiac phase (e.g., when the heart is at rest) is accumulatedfrom multiple source rotations and then LT reconstruction is performedfrom that data. Moreover, the problem of 4D CT imaging can beapproximately solved. Using again the ECG data, the complete cardiaccycle is split into several segments, accumulates the cone beam data foreach segment from multiple source rotations, and then applies LTseparately for each segment. See F. Koo, H. Kudo and L. Zeng,Proceedings of the VIIIth International Conference on FullyThree-Dimensional Reconstruction in Radiology and Nuclear Medicine, SaltLake City, Utah, University of Utah, (2005), where similar techniquesare used together with more conventional inversion algorithms.

Because of the flexibility of LT, its relative stability with respect toinconsistencies in the data (compared with “global” algorithms) and theability to accurately reconstruct edges inside objects, it is expectedthat LT can become a valuable tool, which provides important informationcomplementing well-established inversion techniques.

While the invention has been described, disclosed, illustrated and shownin various terms of certain embodiments or modifications which it haspresumed in practice, the scope of the invention is not intended to be,nor should it be deemed to be, limited thereby and such othermodifications or embodiments as may be suggested by the teachings hereinare particularly reserved especially as they fall within the breadth andscope of the claims here appended.

1. A method of reconstructing an image from cone beam (CB) projectiondata provided by at least one detector, comprising the steps of:collecting CB projection data of an object; storing the collected CBprojection data in a memory; and reconstructing the image of the objectfrom a local CB projection data including the steps of: filtering thelocal cone beam projection data using only derivative operations; andbackprojection updating of the image to reconstruct the image withartifacts having a strength that is weaker than the strength of theuseful features, the useful features being high spatial contrastfeatures in the image which reflect the corresponding high spatialcontrast features in the object.
 2. The method of claim 1, wherein thereconstructing step comprises the step of: finding a combination of onlyderivatives of the CB projection data that weaken the artifacts.
 3. Themethod of claim 2, wherein the finding a combination of derivatives ofthe CB projection data step comprises the steps of: collecting cone beamdata that represents a collection of integrals with or without weight ofan unknown function that represents the object being scanned, whereinsaid integrals are over lines and said lines intersect a curve; andfinding the combination of derivatives of the CB projection data that isequivalent to differentiating the data along a direction which istangent to a source trajectory.
 4. The method of claim 3, furthercomprising the step of: computing the combination of derivatives of theCB projection data, which is equivalent to differentiating the dataalong the direction of the tangent of the source trajectory.
 5. Themethod of claim 4, wherein the differentiation step comprises the stepof: computing derivatives${\frac{\partial^{2}}{\partial\alpha^{2}}{D_{f}\left( {s,\left( {\alpha,w} \right)} \right)}},{\frac{\partial^{2}}{{\partial\alpha}{\partial w}}{D_{f}\left( {s,\left( {\alpha,w} \right)} \right)}},{\frac{\partial^{2}}{\partial w^{2}}{D_{f}\left( {s,\left( {\alpha,w} \right)} \right)}}$each CB projection data, wherein α is the angular parameter on thedetector, w is the vertical parameter on the detector, s is theparameter along the source trajectory, x is the point where imagereconstruction is performed, D_(ƒ)(s,(α, w)) is the CB projection data,and ƒ is the function representing the object being scanned.
 6. Themethod of claim 5, wherein the reconstruction step further comprises thesteps of: multiplying the derivatives by the quantities cos² α, −2 cosα(w sin α−h/(2π)), and (w sin α−h/(2π))², respectively; adding themultiplied derivatives to produce a combination of derivatives;multiplying the combination of derivatives by φ(α,w)cos² a to producethe function Φ(s; α, w) according to${{\Phi\left( {{s;\alpha},w} \right)} = {{\varphi\left( {\alpha,w} \right)}\cos^{2}{\alpha\left\lbrack {{\cos^{2}\alpha\frac{\partial^{2}}{\partial\alpha^{2}}} - {2\cos\;{\alpha\left( {{w\;\sin\;\alpha} - {h/\left( {2\pi} \right)}} \right)}\frac{\partial^{2}}{{\partial\alpha}{\partial w}}} + {\left( {{w\;\sin\;\alpha} - {h/\left( {2\pi} \right)}} \right)^{2}\frac{\partial^{2}}{\partial w^{2}}}} \right\rbrack}{D_{f}\left( {s,\left( {\alpha,w} \right)} \right)}}},$computing the CB local tomography function${{g(x)} = {\int_{I}{\frac{R^{2}}{V^{2}\left( {s,x} \right)}{\Phi\left( {{s;{\alpha\left( {s,x} \right)}},{w\left( {s,x} \right)}} \right)}{{\mathbb{d}s}.}}}},$wherein the source trajectory is a helix, and R is the radius of thehelix, h is the pitch of the helix, α(s, x), w(s, x) are, respectively,the α and w coordinates of the point x projected onto detector, whichcorresponds to the point s of the source trajectory, φ(s, x)is a cut-offfunction which determines the portion of the source trajectory that isused for image reconstruction at x, φ(α, w) is the function determinedby the relation φ(α(s, x), w(s, x))=φ(s, x),${{V\left( {s,x} \right)} = {R - \frac{{x_{1}{y_{1}(s)}} + {x_{2}{y_{2}(s)}}}{R}}},$x₁, x₂ are two coordinates of the point x, and y₁(s), y₂(s) are twocoordinates of the point s of the source trajectory.
 7. The method ofclaim 4, further comprising the step of: processing the derivativeresult by multiplying it by a weight to produce a processed data.
 8. Themethod of claim 1, further comprising the step of: determining a shiftbetween a location of the reconstruction point shown on thereconstructed image and an actual location of a corresponding usefulfeature of the object to determine a location error of a moving objectto correct for the location of the useful features of the object.
 9. Asystem for reconstructing an image of an object comprising: a scannerfor scanning the object in three-dimensions to produce a cone beamprojection data; a processing unit having a memory for storing the conebeam projection data and executing instructions; a first set ofinstructions for calculating derivatives of the cone beam projectiondata; a second set of instruction for reconstructing useful features ofthe image from a local cone beam projection data including: a firstsubset of instructions to filter the local cone beam projection dataconsisting only of derivative operations; and a second subset ofinstructions to backprojection update the filtered local cone beamprojection data to reconstruct the image with suppressed artifacts thatare weaker than the strength of the useful features, the useful featuresbeing high spatial contrast features in the reconstructed imagecorresponding to the high spatial contrast features in the object; and adisplay connected with said processing unit for displaying thereconstructed image.
 10. The system of claim 9, wherein the first subsetof instructions includes instructions for differentiating the cone beamprojection data along a direction tangent to the curve that representsthe scanner trajectory.
 11. The system of claim 9, wherein the first setof instructions comprises: a third subset of instructions for finding adirection of the derivative of the cone beam projection data that willresult in suppressed artifacts; and a fourth subset of instructions fordifferentiating the cone beam projection data according to thedirection.
 12. The system of claim 9, further comprising: a third set ofinstructions for determining a shift between a location of thereconstruction point shown on the reconstructed image and an actuallocation of the point on the object to correct for the location of theuseful features of the actual image.
 13. A method for cone beam localtomography comprising the steps of: loading current cone beam projectiondata of an object into a computer memory, wherein the cone beamprojection data corresponds to a focus of beams of radiation located aty(s), wherein y(s) is a curve along which a source of radiation movesrelative to an object being scanned; calculating derivatives of the conebeam projection data in a direction that results in suppression of theartifacts; processing the derivative results of a local cone beamprojection data to reconstruct the image with suppressed artifacts thathave a strength weaker than the strength of the useful features that arehigh spatial contrast features corresponding to the high spatialcontrast features of the object; applying back projection to theprocessed data; and displaying the reconstructed image with suppressedartifacts without suppressing the strength of useful features.
 14. Themethod of claim 13, wherein the differentiation step comprises the stepsof: determining a direction of the derivative that results insuppression of the artifacts, wherein the direction is tangent to thecurve y(s); and calculating derivatives of the cone beam projection dataalong the direction tangent to the curve y(s) to suppress artifacts. 15.A method of reconstructing an image from local cone beam (CB) dataprovided by at least one detector, comprising the steps of: collectingCB projection data of an object; storing the collected CB projectiondata in a memory; reconstructing the image from the local CB projectiondata, wherein for each reconstruction point x the local cone beam isdefined as data corresponding to rays of radiation passing through asmall neighborhood of the point x, the size of the small neighborhoodnot exceeding that which is required to numerically compute a derivativeat each reconstruction point; and wherein artifacts of the reconstructedimage have a weaker strength than the strength of the useful features,the useful features being high spatial contrast features in thereconstructed image that correspond to the high spatial contrastfeatures in the object.